      SUBROUTINE E01DAZ(NXORDR,NXKNTS,XMIN,XMAX,NYORDR,NYKNTS,YMIN,YMAX,
     *                  MX,X,MY,Y,F,IF1,IF2,LF,XLAM,LXLAM,YLAM,LYLAM,C,
     *                  IC1,IC2,LC,XUFCTR,LXU,YUFCTR,LYU,E,IE1,IE2,LE,
     *                  XKNOT,LXKNOT,YKNOT,LYKNOT,XROW,LXROW,YROW,LYROW,
     *                  WRK,LWRK,IFAIL)
C     MARK 14 RELEASE. NAG COPYRIGHT 1989.
C     A DASL routine, unmodified except for name.
C
C     **********************************************************
C
C     D A S L  -  DATA APPROXIMATION SUBROUTINE LIBRARY
C
C     SUBROUTINE B2IR      BIVARIATE POLYNOMIAL SPLINE
C     ===============      INTERPOLANT TO DATA ON
C                          RECTANGULAR MESH
C                          - BASE VERSION
C
C     CREATED 24 09 81.  UPDATED 15 11 82.  RELEASE 00/07
C
C     AUTHOR ... MAURICE G. COX.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND.
C
C     (C)  CROWN COPYRIGHT 1981-1982
C
C     **********************************************************
C
C     B2IR.  DETERMINES A POLYNOMIAL SPLINE
C     INTERPOLANT TO DATA SPECIFIED AT ALL MESHPOINTS OF A
C     RECTANGULAR MESH.
C
C     INPUT PARAMETERS
C        NXORDR   ORDER (DEGREE + 1) OF  SPLINE  S  IN  X
C        NXKNTS   NUMBER OF INTERIOR X-KNOTS
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF X-INTERVAL
C        NYORDR   ORDER OF  S  IN  Y
C        NYKNTS   NUMBER OF INTERIOR Y-KNOTS
C        YMIN,
C        YMAX     LOWER AND UPPER ENDPOINTS OF Y-INTERVAL
C        MX       NUMBER OF X-MESHLINES
C        X        X-MESHLINE VALUES
C        MY       NUMBER OF Y-MESHLINES
C        Y        Y-MESHLINE VALUES
C        F        DEPENDENT VARIABLE VALUES
C        IF1,
C        IF2      INDEX INCREMENTS OF  F
C        LF       DIMENSION OF  F
C
C     OUTPUT (AND ASSOCIATED DIMENSION) PARAMETERS
C        XLAM     INTERIOR X-KNOTLINES
C        LXLAM    DIMENSION OF  XLAM.  .GE. MAX(NXKNTS, 1).
C        YLAM     INTERIOR Y-KNOTLINES
C        LYLAM    DIMENSION OF  YLAM.  .GE. MAX(NYKNTS, 1).
C        C        B-SPLINE COEFFICIENTS OF  S
C        IC1,
C        IC2      INDEX INCREMENTS OF  C
C        LC       DIMENSION OF  C
C
C     WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C        XUFCTR   BAND UPPER TRIANGULAR MATRIX FOR  X
C        LXU      DIMENSION OF  XUFCTR.
C                    .GE. NXORDR*(MX + NXKNTS + 1)/2
C        YUFCTR   BAND UPPER TRIANGULAR MATRIX FOR  Y
C        LYU      DIMENSION OF  YUFCTR.
C                    .GE. NYORDR*(MY + NYKNTS + 1)/2
C        E        INTERMEDIATE COEFFICIENTS FROM FIRST FIT
C        IE1,
C        IE2      INDEX INCREMENTS OF  E
C        LE       DIMENSION OF  E  .GE.  MY*MX.
C        XKNOT    KNOTS USED IN COMPUTING B-SPLINE BASIS
C                    FOR CURRENT X-VALUE
C        LXKNOT   DIMENSION OF  XKNOT.  .GE. 2*NXORDR.
C        YKNOT    KNOTS USED IN COMPUTING B-SPLINE BASIS
C                    FOR CURRENT Y-VALUE
C        LYKNOT   DIMENSION OF  YKNOT.  .GE. 2*NYORDR.
C        XROW     OBSERVATIONAL ROW CONTAINING B-SPLINE
C                    BASIS VALUES FOR CURRENT X-VALUE
C        LXROW    DIMENSION OF  XROW.  .GE. MX.
C        YROW     OBSERVATIONAL ROW CONTAINING B-SPLINE
C                    BASIS VALUES FOR CURRENT Y-VALUE
C        LYROW    DIMENSION OF  YROW.  .GE. MY.
C        WRK      WORKSPACE REQUIRED BY  ERB
C        LWRK     DIMENSION OF  WRK.  .GE. MAX(MX, MY).
C
C     FAILURE INDICATOR PARAMETER
C        IFAIL    FAILURE INDICATOR
C                    1 - NUMERICALLY SINGULAR SYSTEM
C
C     ----------------------------------------------------------
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN, YMAX, YMIN
      INTEGER           IC1, IC2, IE1, IE2, IF1, IF2, IFAIL, LC, LE, LF,
     *                  LWRK, LXKNOT, LXLAM, LXROW, LXU, LYKNOT, LYLAM,
     *                  LYROW, LYU, MX, MY, NXKNTS, NXORDR, NYKNTS,
     *                  NYORDR
C     .. Array Arguments ..
      DOUBLE PRECISION  C(LC), E(LE), F(LF), WRK(LWRK), X(MX),
     *                  XKNOT(LXKNOT), XLAM(LXLAM), XROW(LXROW),
     *                  XUFCTR(LXU), Y(MY), YKNOT(LYKNOT), YLAM(LYLAM),
     *                  YROW(LYROW), YUFCTR(LYU)
C     .. Local Scalars ..
      INTEGER           IERROR
C     .. External Subroutines ..
      EXTERNAL          E01DAW
C     .. Executable Statements ..
C
C     CONSTRUCT POLYNOMIAL SPLINE INTERPOLANT TO F-VALUES
C     ALONG EACH Y-MESHLINE
C
      CALL E01DAW(NXORDR,XMIN,XMAX,MX,X,MY,F,IF1,IF2,LF,NXKNTS,XLAM,
     *            LXLAM,E,IE1,IE2,LE,XUFCTR,LXU,XKNOT,LXKNOT,XROW,LXROW,
     *            WRK,LWRK,IERROR)
      IF (IERROR.NE.0) GO TO 20
C
C     CONSTRUCT POLYNOMIAL SPLINE INTERPOLANT TO THE
C     COEFFICIENTS SO OBTAINED ALONG EACH X-MESHLINE
C
      CALL E01DAW(NYORDR,YMIN,YMAX,MY,Y,MX,E,IE2,IE1,LE,NYKNTS,YLAM,
     *            LYLAM,C,IC2,IC1,LC,YUFCTR,LYU,YKNOT,LYKNOT,YROW,LYROW,
     *            WRK,LWRK,IERROR)
   20 IFAIL = IERROR
      RETURN
C
C     END E01DAZ
C
      END
